Velocity Profiles in Repulsive Athermal Systems under Shear 
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We conduct molecular dynamics simulations of athermal systems undergoing boundary-driven 
planar shear flow in two and three spatial dimensions. We find that these systems possess nonlinear 
mean velocity profiles when the velocity u of the shearing wall exceeds a critical value u c . Above 
u c , we also show that the packing fraction and mean-square velocity profiles become spatially- 
dependent with dilation and enhanced velocity fluctuations near the moving boundary. In systems 
with overdamped dynamics, u c is only weakly dependent on packing fraction <j>. However, in systems 
with underdamped dynamics, u c is set by the speed of shear waves in the material and tends to 
zero as (f> approaches <f) c . In the small damping limit, <f) c approaches values for random close-packing 
obtained in systems at zero temperature. For underdamped systems with (j> < <f) c , u c is zero and 
thus they possess nonlinear velocity profiles at any nonzero boundary velocity. 

PACS numbers: 83.10.Rs,83.50.Ax,45.70.Mg,64.70.Pf 



Driven, dissipative systems are ubiquitous in nature 
(occurring much more frequently than equilibrium ther- 
mal systems) and display complex behaviors such as hys- 
tcretic and spatially-dependent flows. Moreover, there 
is no complete theoretical description for these systems, 
which makes it difficult to predict how they will respond 
to applied loads such as shear stress. Many of these sys- 
tems such as granular materials Q, Q, metallic glasses 
and complex fluids, for example emulsions Q, foams 
[a, [fj, and worm- like micelles 0, do not flow homo- 
geneously with a linear velocity profile when they are 
sheared. Shear localization or banding can occur where 
a small fraction of the system near one of the bound- 
aries undergoes strong shear flow while the remainder of 
the system is nearly static. Despite much intense work, 
a complete description of how these systems respond 
to shear stress is not available. We perform molecular 
dynamics simulations of repulsive athermal particulate 
systems in two (2D) and three (3D) spatial dimensions 
undergoing boundary-driven planar shear flow to study 
mechanisms that give rise to spatially inhomogeneous ve- 
locity profiles. 

We will answer several important questions in this 
letter. First, does the packing fraction of the system 
strongly influence the shape of the velocity profiles? 
Most previous simulations investigating velocity profiles 
in sheared systems have been performed either near ran- 
dom close-packing as in simulations of granular materials 
Q or at high density as in studies of Lennard- Jones liq- 
uids and glasses [Tol ITU . However, a systematic study 
of the role of density has not been performed. Nonlinear 
mean velocity profiles have been found at both high den- 
sity and near random close packing, but it is not clear 
whether the same physical mechanism is responsible in 
both regimes. 

We also consider the influence of the speed u of the 
shearing boundary on the mean velocity profiles. Re- 
sults from previous simulations of glassy systems [HI Hill 
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FIG. I: Shear stress E vs. velocity u of the wall moving at 
constant velocity (circles) and constant stress (squares) for 
a 2D underdamped system with linear spring interactions at 
4> = 0.85. Ti vs is the yield stress at constant stress. Previous 
studies [Tol ITU have shown that mean velocity profiles can 
switch from nonlinear to linear as u increases above no, where 



uo is the wall velocity at E 



We will show below 



that the mean velocity profiles become nonlinear again when 
u > u c , where u c > uq is always satisfied. 

indicate that a critical velocity uq exists below which the 
mean velocity profiles become nonlinear [T2| . In these 
systems, the yield stress £j, s at constant shear stress 
(shear stress required to move the boundary at nonzero 
velocity) is larger than the yield stress Y, yv at constant 
velocity (shear stress in the limit of zero velocity) . When 
the shear stress is between these two values, part of the 
system can flow while other parts remain nearly static. 
In this letter, we concentrate instead on the larger u 
regime, and ask whether the velocity profiles remain lin- 
ear for all u > uq ■ We will show that another transition 
takes place — the velocity profile switches from linear to 
nonlinear — when the boundary velocity exceeds u c > uq. 
The onset of nonlinear velocity profiles at large u coin- 
cides with the appearance of nonuniform packing fraction 
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FIG. 2: (a) Average velocity (v x ) (normalized by u) in the flow 
direction, (b) local packing fraction <j>, and (c) velocity fluctu- 
ations (5vx, y ) in the x- (solid lines) and y-directions (symbols) 
as a function of height y/L y from the stationary wall in a 2D 
system with harmonic spring interactions and underdamped 
dynamics (b* = 0.01) at tj> = 0.85. In each panel, 4 boundary 
velocities are shown; triangles, diamonds, squares, and circles 
correspond to u = 0.075, 0.15, 0.37, and 0.75, respectively. 
The inset to (c) compares velocity fluctuations in the x- and 
y- direct ions at u — 0.75. 

and temperature profiles. Both u regimes are depicted in 
Fig. fusing the flow curve for an underdamped athermal 
system in 2D. 

In order to demonstrate these results, we performed a 
series of molecular dynamics simulations of soft repulsive 
athermal systems undergoing boundary-driven shear flow 
under conditions of fixed volume, number of particles N, 
and velocity of the top shearing wall u. The systems were 
composed of N/2 large particles and N/2 small particles 
with equal mass m and diameter ratio 1.4 to prevent crys- 
tallization and segregation. Initial states were prepared 
by quenching the system from random initial positions 
to zero temperature ^3| using the conjugate gradient 
method 0] to minimize the system's total potential en- 
ergy. During the quench, periodic boundary conditions 
were implemented in all directions. Following the quench, 
particles with y-coordinates y > L y (y < 0) were chosen 
to comprise the top (bottom) boundary. The walls were 



therefore rough and amorphous. Results did not depend 
on the thermal quench rate provided the systems were 
sheared long enough to remove initial transients. 

Shear flow in the x-direction with a shear gradient in 
the y-direction was created by moving all particles in the 
top wall at fixed velocity u in the x-direction relative to 
the stationary bottom wall. Therefore, particles in the 
walls do not possess velocity fluctuations. During the 
shear flow, periodic boundary conditions were imposed 
in the x- and z-directions (in 3D). The system-size was 
varied in the range N = [256,3072] to assess finite-size 
effects. Only small sample sizes were required in the x 
and z directions. In contrast, more than « 50 particle 
layers were required in the shear-gradient direction to 
remove finite-size effects. Most simulations were carried 
out using L x — L z — 18a and L y = 72a, where a is the 
small particle diameter. The systems were sheared for a 
strain of 5 to remove initial transients and then quantities 
like velocity, pressure and shear stress (obtained from 
the microscopic pressure tensor 15J), and local packing 
fraction were measured as a function of distance y from 
the stationary wall. Averaged quantities were obtained 
by sampling between strains of 5 to 10. 

Bulk and boundary particles interact via the follow- 
ing pairwise, finite-range, purely repulsive potential: 
V(rij) — e (1 — rij/aij) a /a, where a = 2, 5/2 cor- 
respond to harmonic and Hertzian spring interactions, 
e is the characteristic energy scale of the interaction, 
c%j = (<7j + aj)/2 is the average diameter of particles 
i and j, and r*y is their separation. The interaction po- 
tential is zero when ry > . Our results were obtained 
over a range of packing fraction from = [0.58, 0.80] in 
3D and <j> — [0.75, 1.0] in 2D, which allows us to probe 
packing fractions both above and below random close- 
packing [l3|. The units of length, energy, and time are 
a, e, and a^/m/e, respectively. 

For athermal or dissipative dynamics, the position and 
velocity of each particle are obtained by solving 0] 



ni ■ 



dt 2 



F r 



Vj), 



(1) 



where F( = — dV{rij) / dr^rij , the sums over j only 
include particles that overlap i, Vi is the velocity of parti- 
cle i, and b > is the damping coefficient. Although for 
the present discussion we neglect the particles' rotational 
degrees of freedom, we have shown that including these 
does not qualitatively change any of our results [17| • The 
dynamics can be changed from underdamped to over- 
damped by tuning the dimensionless damping coefficient 
b* = ba/^/em above b* = \[2. Frictionless granular ma- 
terials and model foams can be studied using b* < b* 0] 
and b* » b* c [3, respectively. 

Three physical parameters, the packing fraction <j), the 
velocity u of the moving boundary, and the dimension- 
less damping coefficient &*, strongly influence the shape 
of the mean velocity profile. First, we find that a criti- 
cal boundary velocity u c exists that separates linear from 
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FIG. 3: Critical velocity u c of the moving wall versus pack- 
ing fraction cj> in (a) 2D and (b) 3D systems with harmonic 
(circles) and Hertzian (squares) spring interactions. The open 
and filled symbols correspond to b* — 0.01 and b* = 5, respec- 
tively. For underdamped systems, we plot ut/2 for harmonic 
(small circles) and Hertzian (small squares) spring interac- 
tions, where ut is the shear wave speed. 

nonlinear flow behavior. For u < u c (but not in the qua- 
sistatic flow regime), the mean velocity profiles in the 
flow direction are linear; however, when u > u c they 
become nonlinear. The width of the shearing region de- 
creases as u continues to increase above u c . This is shown 
in Fig. EJa) for an underdamped (b* <C b*) system in 
2D with harmonic spring interactions at <f> = 0.85. As 
u is increased above u c ps 0.08, the mean velocity pro- 
files (v x (y)} become more and more nonlinear. When 
the boundary velocity has increased to u = 0.75, ap- 
proximately 80% of the system is nearly static, while the 
remaining 20% undergoes shear flow. 

We also monitored the local packing fraction and 
mean-square velocity fluctuations (or kinetic tempera- 
ture) during shear. These are shown for the same dense 
system with underdamped dynamics in Fig. [5] (b) and 
(c). We find that when the mean velocity profile is lin- 
ear, the packing fraction and velocity fluctuations are 
spatially uniform. Moreover, the velocity fluctuations 
in the x- and y-directions are identical. However, when 
the boundary velocity exceeds u c , the packing fraction 
and mean-square velocity profiles become spatially de- 
pendent. In this regime, the compressional forces induced 
by the shearing boundary are large enough to cause dila- 
tancy. The system becomes less dense near the shearing 
wall and more compact in the nearly static region. In 



addition, the shearing wall induces a kinetic tempera- 
ture gradient with velocity fluctuations larger near the 
shearing boundary. The kinetic temperature also be- 
comes anisotropic with (Svl) < (Svy) when u > u c . Thus, 
several phenomena occur simultaneously as the boundary 
velocity is increased above u c : 1) the velocity profile be- 
comes nonlinear, 2) the system dilates near the shearing 
boundary and compacts in the bulk, and 3) the kinetic 
temperature becomes higher near the shearing wall. 

We have measured the critical wall velocity ztc as a 
function of packing fraction (j) for systems with under- 
damped dynamics and harmonic and Hertzian spring in- 
teractions in 2D and 3D. These measurements are shown 
in Fig. |3 We find that u c is nearly constant at large (j> but 
then decreases sharply as <f> approaches a critical packing 
fraction <p c . For <f> < </> c , u c = with c/> c « 0.81 — 0.82 
in 2D for harmonic and Hertzian springs and 4> c f=a 0.61 
in 3D for harmonic springs. We expect that Hertzian 
springs will give a similar result for (j) c in 3D. These val- 
ues for <p c are close to recent measurements of random 
close-packing in systems at zero temperature Il3|. We 
have also measured u in underdamped systems |17|, uq 
decreases strongly near random close packing, but for all 
systems studied uq < u c . In particular, we find that 
when u c = 0, uq = also. 

A possible interpretation of the critical wall velocity u c 
can be obtained by comparing the time it takes the sys- 
tem to shear a unit strain to the time it takes a shear 
wave (with speed ut) to traverse the system and re- 
turn to the shearing boundary. This simple argument 
predicts u c — Ut/2. ut can be obtained by studying 
the transverse current correlation function CV(w, A;) as a 
function of frequency w and wavenumber k — 2irna/L x 
(n =integer) and the resulting dispersion relation o>t(&) 
[19J. In Fig. |21 we compare ut = dtox/dk (for n = 3 to 
12) and u c as a function of <fi for both potentials in 2D 
and for harmonic springs in 3D [2(J . Although deviations 
occur close to <f> c , we find that u c agrees very well with 
ut/2 over a wide range of <f>. 

What is the shape of mean velocity profiles in dilute 
underdamped systems with (j) < <j) c l Since u c is zero 
below (j> c , we expect that mean velocity profiles in these 
dilute systems are nonlinear for all nonzero u. This is 
indeed what we find for all systems studied. Fig.0|shows 
the mean velocity profiles for a 2D underdamped system 
at (f> < <j) c over three decades in u. In contrast to the 
behavior in dense systems, the velocity profiles are not 
monotonic in u. However, there is a range of boundary 
velocities (one decade) over which the velocity profiles 
collapse onto a common exponential profile. A robust 
exponential profile has also been found over a wide range 
of shear rates in experiments of granular materials 2}. 
Similarly to the systems characterized by large <f>, low 
4> ones are accompanied by spatially-dependent packing 
fraction and mean-square velocity profiles. 

Boundary-driven shear flow in overdamped systems 
is, however, substantially different from that in under- 
damped systems since velocities of neighboring particles 
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FIG. 4: Average velocity {v x )/u in the shear flow direction 
as a function of height y/L y from the stationary wall in a 2D 
underdamped system (b* = 0.01) at <j> < cj> c . Three boundary 
velocities are shown; squares, downward triangles, and pluses 
correspond to u — 0.38, 0.038, and 7.7 x 10 -4 , respectively. 
The inset shows that there is a wide range of u from 0.0077 
(leftward triangles) to 0.077 (circles) over which the velocity 
profiles collapse. 

are strongly coupled. Fig. [21 shows that in the over- 
damped limit (b* 3> &*), the critical boundary velocity is 
nearly independent of <fr over the studied range in both 
2D and 3D. We also find that u c increases linearly with 
b* , thus the velocity profiles tend toward linear profiles 
as the dissipation increases at fixed u. 

In this letter we present results of molecular dynam- 
ics simulations of repulsive athermal systems undergoing 
boundary-driven shear flow in 2D and 3D. We demon- 
strate that a critical boundary velocity u c exists (at large 



u) that signals the onset of spatial inhomogeneity. When 
u exceeds u c , the mean velocity profiles become nonlin- 
ear, the system becomes dilated near the moving wall 
and compressed near the stationary wall, and the sys- 
tem possesses a nonuniform kinetic temperature profile 
with higher temperature near the moving wall. For un- 
derdamped systems, u c is nearly constant at large <f> but 
decreases strongly at lower <fi until it vanishes at <fr c . (f> c 
depends on the damping coefficient but approaches re- 
cent estimates of random close-packing [13j in the small 
damping limit. In this limit, u c is determined by the 
speed of shear waves in the material, ut- When u ex- 
ceeds u c = Ut/2, large shear strain occurs before shear 
waves are able to traverse the system. In the overdamped 
limit, u c is nearly independent of <f> over the studied range 
and scales linearly with the damping coefficient b*. 

Possible future directions include adding friction to the 
model in Eq.^ which would allow us to make more direct 
contact with experiments on granular materials P, Q{. 
Are nonlinear velocity profiles more likely to occur in 
systems with frictional particles? How does the propaga- 
tion of shear waves change in that case? We are currently 
investigating these important questions. 
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